#' @title RDD poliynomial local regression
#' @description no parametric regression discontinuity using polynomial local regression estimator with and without sample weights
#' @param D variable with: 1 for treatment group and 0 for control group
#' @param Z continuous variable
#' @param Y outcome variable
#' @param c cutpoint
#' @param type sharp or fuzzy
#' @param TreatLeft T if the treatment is applied for below of the cutpoint
#' @param W sample weights 1/pik, NULL if there aren't sample weights
#' @param svydesign_ design sample using survey library
#' @param vcov_ variances and covariances type: HCO, HC1, HC2 y HC3
#' @param kernel_ kernel type: triangular, epanechnikov, tricube, uniform, gaussian y cosine
#' @param h bandwidth
#' @param q order polynomial
#' @return Impact estimation, e.e. and t test
#' @examples
#'
RDpol <- function(D = muestra1$D, Z = muestra1$x, Y = muestra1$y,c =0,type="sharp",TreatLeft=T,weights=NULL,svydesign_=muestra,vcov_,kernel_="uniforme",h=2,q=2){
library(survey)
if(!is.null(weights)){
muestra <- data.frame(D,Z,Y,weights)
}else{
muestra <- data.frame(D,Z,Y)
}
muestra$intercepto <- 1
muestra$Z_C <- muestra$Z-c
muestra$ZD <- muestra$Z*muestra$D
if(kernel_!="gaussian"){
muestra <- muestra[which(muestra$Z>=c-h & muestra$Z<=c+h),]
if(!is.null(weights)){
svydesign_ <- subset(svydesign_,Z>=c-h & Z<=c+h)
}
}
Y <- muestra$Y
if(!is.null(weights)){
W <- kernel_sRD(kernel_=kernel_,weights=1/muestra$weights,h=h,X=muestra$Z,c=c)
}else{
W <- kernel_sRD(kernel_=kernel_,weights=NULL,h=h,X=muestra$Z,c=c)
}
if(type!="sharp"){
muestra$T_ <- ifelse(muestra$Z<c,1,0)
if(TreatLeft==F){
muestra$T_ <- 1-muestra$T_
}
muestra$ZT_ <- muestra$Z*muestra$T_
}
if(type=="sharp"){
Xm1 <- as.matrix(muestra[,c("intercepto","Z_C","D","ZD")])
for(i in 2:q){
Xm1 <- cbind(Xm1,Xm1[,2]^i,Xm1[,4]^i)
colnames(Xm1)[(3+2*(i-1)):(4+(2*(i-1)))] <- c(paste0("Z_C_",q),paste0("ZD_",q))
}
impacto.1_2 <- solve(t(Xm1)%*%W%*%Xm1)%*%(t(Xm1)%*%W%*%Y)
beta <- impacto.1_2[3,]
}else{if(type=="fuzzy"){
Zm1 <- as.matrix(muestra[,c("intercepto","Z_C","T_","ZT_")])
for(i in 2:q){
Zm1 <- cbind(Zm1,Zm1[,2]^i,Zm1[,4]^i)
colnames(Zm1)[(3+2*(i-1)):(4+(2*(i-1)))] <- c(paste0("Z_C_",q),paste0("ZD_",q))
}
Xm1 <- as.matrix(muestra[,c("intercepto","Z_C","D","ZT_")])
for(i in 2:q){
Xm1 <- cbind(Xm1,Xm1[,2]^i,Xm1[,4]^i)
colnames(Xm1)[(3+2*(i-1)):(4+(2*(i-1)))] <- c(paste0("Z_C_",q),paste0("ZD_",q))
}
beta1stg <- solve(t(Zm1)%*%W%*%Zm1)%*%(t(Zm1)%*%W%*%muestra$D)
muestra$Xest <- Zm1%*%beta1stg
Xest <- as.matrix(muestra[,c("intercepto","Z_C","Xest","ZT_")])
for(i in 2:q){
Xest <- cbind(Xest,Xest[,2]^i,Xest[,4]^i)
colnames(Xest)[(3+2*(i-1)):(4+(2*(i-1)))] <- c(paste0("Z_C_",q),paste0("ZD_",q))
}
impacto.1_2 <- solve(t(Xest)%*%W%*%Xest)%*%(t(Xest)%*%W%*%Y)
beta <- impacto.1_2[3,]
}
}
# Varianza sin pesos
Yest <- Xm1%*%impacto.1_2
sigma2 <- sum((Y - Yest)^2)/(dim(muestra)[1]-length(impacto.1_2))
M_ <- (Y - Yest)^2
M <- matrix(0,length(M_),length(M_))
diag(M) <- M_
if(type=="sharp"){
if(is.null(weights)){
if(vcov_=="HC0"){
var <- solve(t(Xm1)%*%W%*%Xm1)%*%t(Xm1)%*%W%*%M%*%t(W)%*%Xm1%*%solve(t(Xm1)%*%W%*%Xm1)
}else{if(vcov_=="HC1"){
mat <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])
diag(mat) <- (dim(Xm1)[1]/(dim(Xm1)[1]-length(impacto.1_2)))
var <- solve(t(Xm1)%*%W%*%Xm1)%*%t(Xm1)%*%W%*%mat%*%M%*%t(W)%*%Xm1%*%solve(t(Xm1)%*%W%*%Xm1)
}else{if(vcov_=="HC2"){
mat <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])
diag(mat) <- 1/((1-diag(Xm1%*%solve(t(Xm1)%*%Xm1)%*%t(Xm1))))
var <- solve(t(Xm1)%*%W%*%Xm1)%*%t(Xm1)%*%W%*%mat%*%M%*%t(W)%*%Xm1%*%solve(t(Xm1)%*%W%*%Xm1)
}else{if(vcov_=="HC3"){
mat <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])
diag(mat) <- 1/((1-diag(Xm1%*%solve(t(Xm1)%*%Xm1)%*%t(Xm1)))^2)
var <- solve(t(Xm1)%*%W%*%Xm1)%*%t(Xm1)%*%W%*%mat%*%M%*%t(W)%*%Xm1%*%solve(t(Xm1)%*%W%*%Xm1)
}}}}
}else{
var1 <- solve(t(Xm1)%*%W%*%Xm1) # revisar si dejo o pongo sigma2
muestra$ek2 <- muestra$Y-Xm1%*%impacto.1_2
diag(M) <- muestra$ek2
mat <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])
if(vcov_=="HC0"){
diag(mat) <- 1
}else{if(vcov_=="HC1"){
diag(mat) <- (dim(Xm1)[1]/(dim(Xm1)[1]-length(impacto.1_2)))
}}
W2 <- kernel_sRD(kernel_=kernel_,weights=NULL,h=h,X=muestra$Z,c=c)
svydesign_$variables$uk <- mat%*%M%*%t(W2)%*%Xm1
varmuestra <- vcov(svytotal(~uk,svydesign_))
var <- (var1)%*%varmuestra%*%(var1)
}
}else{if(type=="fuzzy"){
if(is.null(weights)){
if(vcov_=="HC0"){
var <- solve(t(Xest)%*%W%*%Xest)%*%t(Xest)%*%W%*%M%*%t(W)%*%Xest%*%solve(t(Xest)%*%W%*%Xest)
}else{if(vcov_=="HC1"){
mat <- matrix(0,dim(Xest)[1],dim(Xest)[1])
diag(mat) <- (dim(Xest)[1]/(dim(Xest)[1]-length(impacto.1_2)))
var <- solve(t(Xest)%*%W%*%Xest)%*%t(Xest)%*%W%*%mat%*%M%*%t(W)%*%Xest%*%solve(t(Xest)%*%W%*%Xest)
#var <- (dim(Xest)[1]/(dim(Xest)[1]-2))*var
}else{if(vcov_=="HC2"){
mat <- matrix(0,dim(Xest)[1],dim(Xest)[1])
diag(mat) <- 1/((1-diag(Xest%*%solve(t(Xest)%*%Xest)%*%t(Xest))))
var <- solve(t(Xest)%*%Xest)%*%t(Xest)%*%mat%*%M%*%Xest%*%solve(t(Xest)%*%Xest)
}else{if(vcov_=="HC3"){
mat <- matrix(0,dim(Xest)[1],dim(Xest)[1])
diag(mat) <- 1/((1-diag(Xest%*%solve(t(Xest)%*%Xm1)%*%t(Xest)))^2)
var <- solve(t(Xest)%*%Xest)%*%t(Xest)%*%mat%*%M%*%Xest%*%solve(t(Xest)%*%Xest)
}}}}
}else{
var1 <- solve(t(Xest)%*%W%*%Xest) # revisar si dejo o pongo sigma2
muestra$ek2 <- muestra$Y-Xm1%*%impacto.1_2
diag(M) <- muestra$ek2
mat <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])
if(vcov_=="HC0"){
diag(mat) <- 1
}else{if(vcov_=="HC1"){
diag(mat) <- (dim(Xm1)[1]/(dim(Xm1)[1]-length(impacto.1_2)))
}}
W2 <- kernel_sRD(kernel_=kernel_,weights=NULL,h=h,X=muestra$Z,c=c)
svydesign_$variables$uk <- mat%*%M%*%t(W2)%*%Xm1
varmuestra <- vcov(svytotal(~uk,svydesign_))
var <- (var1)%*%varmuestra%*%(var1)
}
}}
seBeta <- sqrt(var[3,3])
significativo <- ifelse(abs(beta[i]/(seBeta[i]))>=qt(0.975,dim(muestra)[1]-dim(Xm1)[2]),"Significativo","No Significativo")
impact <- data.frame(round(beta,3),round(seBeta,5),pvalue=(1-pt(abs(beta/(seBeta)),df = dim(muestra)[1]-length(impacto.1_2)))*2)
colnames(impact) <- c("Impacto","S.E.","P-Value")
impact$significativo <- ifelse(impact$`P-Value`<0.01,"***",ifelse(impact$`P-Value`<0.05,"**",ifelse(impact$`P-Value`<0.1,"*","")))
colnames(impact)[4] <- ""
print(impact)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.